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■ Abstract 

(N : 

The aim of this paper is to present a kinetic numerical scheme for the computations 
of transient pressurised flows in closed water pipes. Firstly, we detail the mathematical 
model written as a conservative hyperbolic partial differentiel system of equations, and then 
we recall how to obtain the corresponding kinetic formulation. Then we build the kinetic 
scheme ensuring an upwinding of the source term due to the topography performed in a 
close manner described by Perthame et al. [131 [T] using an energetic balance at microscopic 
level. The validation is lastly performed in the case of a water hammer in an uniform pipe: 
• we compare the numerical results provided by an industrial code used at EDF-CIH (France), 

which solves the Allievi equation (the commonly used equation for pressurised flows in pipes) 
by the method of characteristics, with those of the kinetic scheme. It appears that they are 
in a very good agreement. 
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^ ! 1 Introduction 

> , 

The work presented in this article is the first step in a more general project: the modelisation 
of unsteady mixed water flows in open channels and in pipes, its kinetic formulation and its 
■ numerical resolution by a kinetic scheme. 

. Since we are interested in flows occuring in closed pipes, it may happen that some parts of 

the flow are free-surface (this means that only a part of the cross-section of the pipe is filled) 
and other parts are pressurised (this means that all the cross-section of the pipe is filled). Let us 
thus recall the current and previous works about mixed flows in closed water pipes by a partial 
state of the art review. Cunge and Wegner [8] studied the pressurised flow in a pipe as if it were 
a free-surface flow by assuming a narrow slot to exist in the upper part of the tunnel, the width 
of the slot being calculated to provide the correct sonic speed. This approach has been credited 
to Preissmann. Later, Cunge [7] conducted a study of translation waves in a power canal 
containing a series of transitions, including a siphon. Pseudoviscosity methods were employed 
to describe the movement of bores in open-channel reaches. Wiggert [16] studied the transient 
flow phenomena and his analytical considerations included open-channel surge equations that 
were solved by the method of characteristics. He subjected it to subcritical flow conditions. His 
solution resulted from applying a similarity between the movement of a hydraulic bore and an 
interface (that is, a surge front wave). Following Wiggert's model, Song, Cardie and Leung [2] 
developed two mathematical models of unsteady free-surface/pressurised flows using the method 
of characteristics (specified time and space) to compute flow conditions in two flow zones. They 
showed that the pressurised phenomenon is a dynamic shock requiring a full dynamic treatment 
even if inflows and other boundary conditions change very slowly. However the Song models 
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do not include the bore presence in the free-surface zone. Hamam and McCorquodale [TT] 
proposed a rigid water column approach to model the mixed flow pressure transients. This 
model assumes a hypothetical stationary bubble across compression and expansion processes. 
Li and McCorquodale [12] extended the rigid water column approach to allow for the transport 
of the trapped air bubble. Recently Fuamba [9] proposed a model for the transition from a 
free surface flow to a pressurised. He wrote the conservation of mass, momentum and energy 
through the transition point and proposed a laboratory validation of his model. In the last few 
years, numerical models mainly based on the Preissmann slot technique have been developed to 
handle the flow transition in sewer systems. Implementing the Preissmann slot technique has 
the advantage of using only one flow type (free-surface flow) throughout the whole pipe and of 
being able to easily quantify the pressure head when pipes pressurise. Let us specially mention 
the work of Garcia-Navarro, Alcrudo and Priestley [TU] in which an implicit method based on 
the characteristics has been proposed. 

The Saint Venant equations, which are written in a conservative form, are usually used 
to describe free surface flows of water in open channels. As said before, they are also used 
in the context of mixed flows (i.e. either free surface or pressurized) using the artifice of the 
Preissmann slot [15],[6]. On the other hand, the commonly used model to describe pressurized 
flows in pipes is the system of the Allievi equations [15]. This system of 1st order partial 
differential equations cannot be written under a conservative form since this model is derived 
by neglecting some acceleration terms. This non conservative formulation is not appropriate for 
a kinetic interpretation of the transition between the two types of flows since we are not able 
to write conservations of appropriate quantities such as momentum and energy. 

Then, it appears that a conservative model and a kinetic interpretation of it which describes 
pressurised flows in closed water pipes could be of great interest. 

The model used in this article to describe pressurised flows in closed water pipe is very closed 
to the Shallow Water equations, and has been established by the authors in [3]. A second order 
well-balanced finite volume scheme was therein presented. We will recall in section [2] the main 
features of this previous work. 

Another approach for the numerical resolution of Shallow Water equations is to use a kinetic 
formulation [13} [T] . The corresponding scheme appears to have interesting theoretical proper- 
ties: the scheme preserves the still water steady state and involves a conservative in-cell entropy 
inequality. Moreover, this type of numerical approximation leads to an easy implementation. 
The present modelisation of pressurised flows is formally very close to the Shallow Water equa- 
tions and it may be very interesting to propose a kinetic formulation and thus to construct a 
kinetic scheme. 

The model for the unsteady mixed water flows in closed water pipes and a finite volume 
discretisation has been previoulsy studied by the authors [3] and a kinetic formulation has been 
proposed in [5] . We will recall in section [2] the main results and the properties of this kinetic 
formulation that will be useful to show the properties of the numerical kinetic scheme such as 
the preservation of the steady state water at rest, and the positivity of the wetted area. 

Section [3] is devoted to the construction of the kinetic scheme. The upwinding of the source 
term due to the topography is performed in a close manner described by Perthame et al. [13] 
using an energetic balance at microscopic level for the Shallow Water equations. 

Finally, we present in section H] a numerical validation of this study by the comparison 
between the resolution of this model and the resolution of the Allievi equation solved by the 
research code belier used at Center in Hydraulics Engineering of Electricite De France (EDF) 
[17| for the case of critical waterhammer tests. 
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2 The mathematical model and the kinetic formulation 



We derived a conservative model for pressurised flows from the 3D system of compressible Euler 
equations by integration over sections orthogonal to the flow axis. 

2.1 The mathematical model : a "Shallow Water like" system of equations 

The equation for conservation of mass and the first equation for the conservation of momentum 
are: 

d t p + dw( P u) = o (i) 

d t (pu) + div(puU) = F x -d x P (2) 



with the speed vector U = ui + vj + wk = ui + V, where the unit vector i is along the main 
axis, p is the density of the water. We use the Boussinesq linearised pressure law (see jT5]): 

where po is the density at the atmospheric pressure P a and (5 the coefficient of compressibility 
of the water. Exterior strengths F are the gravity g and the friction term Sf which is assumed 
to be given by the Manning-Strickler law (see [15] ) : 

Sf = K u \u\ with K = (3) 

where K s > is the Strickler coefficient, depending on the material, and Rh is the so called 

S 

hydraulic radius given by Rh = — — . S represents the cross-section area of the pipe whereas P m 
is the perimeter of the section. Then Equations (JTJ) — (J23) become: 

d t p + d x (pu) + div^ z) (pV) = 

d t (pu) + d x (pu 2 ) + div( ytZ )(puV) = -pg(d x Z + S /) - . 

Assuming that the pipe is infinitely rigid and has a uniform constant cross-section S, and taking 
averaged values in sections orthogonal to the main flow axis, we get the following system written 
in a conservative form for the unknowns M = p S , D = pSu: 

d t {M) + d x {D) = (4) 

d t (D) + dJ^ + c 2 M) = -gM(d x Z + S f ) (5) 

where c = —== is the speed of sound. A complete derivation of this model, taking into account 

the deformations of the pipe, contracting or expanding sections, and a spatial second order Roe- 
like finite volume method in a linearly implicit version is presented in [3] (see [2] for the first order 
implicit scheme). This system of partial differential equation is formally close to the Shallow 
Water equations where the conservative variables are the wet area and the discharge, thus we 
define an "FS-equivalent" wet area (FS for Free Surface) A and a "FS-equivalent discharge" Q 
through the relations: 

M = pS = po A and D = pSu = poQ . 
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Dividing ©-(0) by po we can write this system under the conservative form: 

8 t U + 8 X F(U) = G(x, U) (6) 

where the unknown state is U = (A,QY, the flux vector is F(U) = (Q, — — + c 2 A) t and the 

source term writes G{x,U) = (0, -gA(d x Z + S/))*. This new set of variables allows a more 
natural treatment of mixed flows (see [3]). 

Let us now recall the main properties of the system ([6]) whose proofs can be found in [5]. 

Theorem 1 The system (Ej) is strictly hyperbolic. It admits a mathematical entropy: 

E(A,Q,Z) = ^ + gAZ + c 2 AlnA (7) 

which satisfies, for smooth enough solution, the entropy inequality : 

d t E + d x [u(E + c 2 A)] <0 . 

Also, for the frictionless pipes (Sf = 0), the system |6p admits a family of smooth steady states 
characterized by the relations: 

Q = Au = C\ , 

v 2 

— + gZ + c 2 lnA = C 2 , 



u 
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where C\ and C 2 are two arbitrary constants. The quantity — + g Z + c 2 \nA is also called the 
total head. 

Remark 1 An easy computation leads to the equality: 

d t E + d x [u{E + c 2 A)] = -uS f = -K\u\ 3 < . 

Thus for a frictionless pipe we obtain an entropy equality whereas the entropy inequality is 
strict as soon as a friction term is considered. 

Let us also remark that the still water steady state for frictionless pipe, namely u = 0, 
satisfies: g Z + c 2 In A = C 2 - 

2.2 The kinetic approach 

We present in this section the kinetic formulation for pressurised flows in closed water pipes 
modelised by the preceding system of partial differential equations (see [3] for more details and 
properties). Let us mention that the following results (namely Theroem [2] and Theorem [3]) are 
only valid for frictionless pipes. Let us consider a smooth real function \ which has the following 
properties: 

X(w) = >0, f X(u)du = 1, f ^xi^duo = 1 . (8) 



We then define the density of particles Ai(t,x,£) by the so-called Gibbs equilibrium: 

M(t,x,0 = ^X^~ U(t,X) 



These definitions allow to obtain a kinetic representation of the system ([6]) by the following 
result (see [5] for the proof). 
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Theorem 2 The couple of functions (A, Q) is a strong solution of the system {SJj if and only 
if Ai satisfies the kinetic equation: 

^M+f.^M-^Z-^M-K^ (9) 

for some collision term K(t,x,£) which satisfies for a.e. (t,x) 

[ Kd£ = , I iKd^ = . 
Jr Jr 

This result is a consequence of the following relations verified by the microscopic equilibrium: 

A = [ M(0<%, (10) 



Q = / £M(£K, (11) 

JR 

91 +c >a = [ eM(o<%. (12) 

A Jr 

This theorem produces a very useful consequence: the nonlinear system ([6]) can be viewed 
as a simple linear equation on a nonlinear quantity Ai for which it is easier to find simple 
numerical schemes with good theoretical properties: it is this feature which will be exploited to 
construct a kinetic scheme. 

Theorem 3 Let A(x, t) > and Q(x, t) be two given functions. 
1. The minimum of the energy: 

£{f) = J (jfiO + c 2 f(0 In(/(0) + gZf(0 + c 2 ln(cv^)/(0) <% , 
under the constraints: 

f>o, I f®dt=A, [ e/m = Q , 

JR JR 

is attained by the function: 

.j, i \ A f £ — u(t, x) 
M(t,x,0 = -x' 



x ( UJ ) = -^=exp(- l 4) . (1.3) 



c \ c 

where x is defined by: 

v(/.A — . 

2. Moreover, the function % defined by \13\) ensures us to have the relation 

£(M) = E(A,Q,Z) 

if A and Q are solution of the pressurised flow equations and the entropy E is defined 
by 0). 

3. The Gibbs equilibrium Ai satisfies the still water steady state equation. 
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3 The kinetic scheme 



The spatial domain is a pipe of length L. The main axis of the pipe is divided in N meshes 
m, =1^-1/2 j x i+i/2^ 1 < * < N , of length hi and center Xj. We denote Ax = mhii<j<7v hi. 
At denotes the time step at time t n and we set t n+ \ = t n + At. 

A? 



The discrete macroscopic unknowns are XJf = I ^ n j with 1 < i < N and < n < n max . 

They represent the mean value of U on the cell rrii at time t n . 

If Z{x) is the function describing the bottom elevation, its piecewise constant representation 
is given by Z(x) = Zi\ m% (x) with Zj = Z{xi) for example. 

Replacing Z by Z and neglecting the collision term K(t,x,t]) in a first step, the Equation 
([9]) in the cell mj writes: 

—M + £ • — X =0 for x G rrii . (14) 
at ox 

This equation is a linear transport equation whose explicit discretisation may be done directly by 
the following way. Denoting for x G rrii, f(t n ,x,£) = 

the maxwellian state associated to A" , and Q", the usual finite volume discretisation of the 
Equation (fT4"|) leads to: 

/r +1 (e) = .M?(e) + ^e as) 

where the fluxes M.. 1 have to take into account the discontinuity of the altitude Z at the cell 
interface Xj+i/2- Indeed, noticing that the fluxes can also be written as: 

the quantity 5A47 L = AiT 1 — , i holds for the discrete contribution of the source term 

1+2 1+2 "~2 

gAd x Z in the system for negative velocities t] < due to the upwinding of the source term. 
Thus 5Ai~ i has to vanish for positive velocity £ > 0, as proposed by the choice of the interface 

fluxes below. 

Let us now detail our choice for the fluxes M. x at the interface. It can be justified by 

using a generalised characteristic method for the Equation (without the collision kernel) 
but we give instead a presentation based on some physical energetic balance. Let us denote 
A~ Z i+ i = Zi + \ — Zi and A + Z i+ i = Zi — Zi + \. In order to take into account the neighboring 
cells by means of a natural interpretation of the microscopic features of the system, we formulate 
a peculiar discretisation for the fluxes in (|15p . computed by the following upwinded formulas: 

reflection 



transmission 



reflection 



M+(0 = M?+i(0%o + A<? + i(-0%<2 ff A+z , I € >o (17) 

~2 l +2" 



transmission 



6 



The effect of the source term is made explicit by treating it as a physical potential. The choices 
()16|) -(fT71) are thus a mathematical formalization to describe the physical microscopic behaviour 
of the system. The contribution of the interface x i+1 / 2 to is given by: 

• the particles in the cell mi at time t n with non negative velocities £ through the term 
A4"(£) I^>o and those of them that are reflected (thus taken into account with velocity 
— £) if their kinetic energy is not large enough to overpass the potential difference i.e. 
£ 2 < 2gA ± Z i ,i: see Figure [2j 

• the particles in the cell mj + i at time t n with a kinetic energy enough to overpass the 
potential difference ( £ 2 > 2gA^Z i+ i ) and speed up or down according to this potential 
jump. It is the transmission phenomenon in classical mechanics as shown in Figure [TJ 

£ <0 



5 >o ^ 


2 -2g AZ 


/ 






Z i + 1 






z . 




1 





X i+l/2 

Figure 1: Transmission 



\ <0 

2 N 




^ Z <2gAZ \ 


Z i + 1 




-\ >o 













X i+1/2 



Figure 2: Reflection 

Since we neglected the collision term, it is clear that / n+1 computed by the discretised kinetic 
equation (|15p is no more a Gibbs equilibrium. Therefore, to recover the macroscopic variables 
A and Q, according to the identities (fT0l) - (fTT1) . we set: 

U ? +1 = ( S?i)=f ( l)f? +1 dt (18) 



Now, we can integrate the discretised kinetic equation (I15p against 1 and £ to obtain the 
macroscopic kinetic scheme: 
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The numerical fluxes are thus defined by the kinetic fluxes as follows: 

± def f f 1 \ ± 



^rt £ ur< m (20) 



Remark 2 

We see immediately that the kinetic scheme (jf9p -(|20 p is wetted area conservative. Indeed, 



let us denote the first component of the discrete fluxes (f20|) (Fa). 



(Fa)%, = [ CMf^(OdC 



An easy computation using the change of variables /i =| £ | 2 — 2gA + Z i+ i in the formulas 
(!16p -(fT7|) defining the kinetic fluxes M.^ 1 allows us to show that: 

• Computing the macroscopic state £/ by the formula (|18p or the fluxes by the formula (|20p 
is not easy if the function % verifying the properties ([8]) is not compactly supported. We 
use instead the function defined by: 

xH = ^=i[_v3,V3]H- ( 21 ) 

A™ 

We get ,M™(£) = % —j= If u n_ cv ^3 U n +CA /3|(C)- Of course, from Theorem^ the property on 

2 c v 3 ' 1 

the microscopic energy and the still water steady state is no more valid but we will prove 
in Theorem H] that our proposed kinetic scheme preserves the still water steady state (for 
the frictionless pipes) and the positivity of the equivalent wetted area. 

• In the case where the friction Sf defined by (J3j) is present, it is added at the macroscopic 
level (11911 as an extra source term. 



We are now able to state the main properties of the kinetic scheme. 



Theorem 4 We choose the function x defined by the formula [21]) and we assume the CFL 
condition 

At max (\uf \ +cV3] < Ax. (22) 
l<i<N V / ~ 

Then 

(i) the kinetic scheme $19\)- (2(^) keeps the pseudo wetted area A™ positive. 

(ii) the kinetic scheme (fiff|)-(f^0| ) preserves the still water steady state, 

uf = 0, gZi + c 2 lnAi = K 



Proof of theorem |4] Let us mention that the CFL condition (|22p is obtained from the linear 
discretised kinetic transport equation (|15p for the particular choice of the function \ defined by 
the formula (I2ip . This condition ensures the positivity of the pseudo wetted area as we show 
below. 
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Since Ai = /" d£, it is sufficient to prove that /" > 0. Writing the microscopic 
Jr 

scheme (fT5|) , (fT6j) , (fT?]) , using the CFL condition ([22]) , and the fact that the function \ that we 
have chosen is compactly supported, one may see that if we suppose that Af > 0, then f™ +1 is 
a sum of non- negative quantities. For the second point, setting u" = 0, we prove easily that in 
the discretised kinetic equation (fT5"j) . we have 

This implies = M? (£), which ensures by definition A^ +1 = A? and = Q™. Thus we 

obtain u™ +1 = 0. 



4 Numerical validation 

We present now numerical results of a water hammer test. The pipe of circular cross-section of 
2 m 2 (the diameter therefater is denoted by 5) and thickness e = 20 cm is 2000 m long. The 
altitude of the upstream end of the pipe is 250 m and the slope is 5°. The Young modulus is 
23 10 9 Pa since the pipe is supposed to be built in concrete. 

The density at the atmospheric pressure po is 1000 kg/m 3 and the coefficient of compressibility 
of the water (3 is 5.0 10~ 10 Pa" 1 . 

The wave speed is thus obtained by the formula (see [T5l formula (2.39)]): 

° = 1086.6 m/s" 1 . (23) 



1 + 



peE 



The total upstream head is 300 m. The initial downstream discharge is 10 m 3 /s and we cut the 
flow in 5 seconds. Let us define the piezometric line by: 

(p (p p ) 

piezo = z + 5 + p with p = . (24) 

Po9 

We present a validation of the proposed scheme by comparing numerical results of the proposed 
model solved by the kinetic scheme with the ones obtained by solving Allievi equations by 
the method of characteristics with the so-called belier code: an industrial code used by the 
engineers of the Center in Hydraulics Engineering of Electricite De France (EDF) [17] . Our 
code is written in Fortran and runs during a few seconds on LinuX, Windows and Macintosh 
operating systems. 

A simulation of the water hammer test was done for a CFL coefficient equal to 0.8 (i.e. 
CFL = 0.8) and a spatial discretisation of 1000 mesh points (the mesh size is equal to 2 m). In 
the Figure El we present a comparison between the results obtained by our kinetic scheme and 
the ones obtained by the "belier" code at the middle of the pipe: the behavior of the piezometric 
line, defined by Equation (I24p . and the discharge at the middle of the pipe. One can observe 
that the results for the proposed model are in very good agreement with the solution of Allievi 
equations. In Figure HI we present the piezometric line for the beginning and the end of the 
simulation. One can see that the peak of pressure (observed in the beginning of the simulation) 
is very well obtained by the kinetic scheme. Thus the strength of the water hammer is very 
good predicted by the proposed numerical scheme. A little smoothing effect, observed at the 
end of the simulation, may be probably due to the first order discretisation type. A second 
order scheme could be implemented naturally and produce a better approximation. 
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5 Conclusion 



As mentionned in the introduction, our goal is to build a kinetic scheme for mixed flows. 
Perthame et al. [TJ [T3] have shown that the kinetic approach is relevant for free surface flows 
in open channels: the resulting kinetic scheme is easily implemented an enjoyed very good 
properties (posivity of the wetted area and discrete entropy inequalities). For pressurised flows, 
we have shown that this approach is also relevant. This allows us to investigate the construction 
of a kinetic scheme for mixed flows in closed water pipes. 
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Figure 3: Comparison between the kinetic scheme and the industrial code belier 
Piezometric line (top) and discharge (bottom) at the middle of the pipe 
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Figure 4: Comparison between the kinetic scheme and the industrial code belier 
ginning of the simulation (top) and end of the simulation (bottom) at the middle of the pipe 
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